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Abstract 

The  gas  dynamics  equations,  coupled  with  a  static  gravitational  field,  admit  the  hydro¬ 
static  balance  where  the  fiux  produced  by  the  pressure  is  exactly  canceled  by  the  gravitational 
source  term.  Many  astrophysical  problems  involve  the  hydrodynamical  evolution  in  a  grav¬ 
itational  field,  therefore  it  is  essential  to  correctly  capture  the  effect  of  gravitational  force  in 
the  simulations.  Improper  treatment  of  the  gravitational  force  can  lead  to  a  solution  which 
either  oscillates  around  the  equilibrium,  or  deviates  from  the  equilibrium  after  a  long  time 
run.  In  this  paper  we  design  high  order  well-balanced  finite  difference  WENO  schemes  to 
this  system,  which  can  preserve  the  hydrostatic  balance  state  exactly  and  at  the  same  time 
can  maintain  genuine  high  order  accuracy.  Numerical  tests  are  performed  to  verify  high 
order  accuracy,  well-balanced  property,  and  good  resolution  for  smooth  and  discontinuous 
solutions. 
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1  Introduction 


In  recent  years,  research  on  well-balanced  numerical  methods  for  the  hyperbolic  equations 
with  source  terms  (also  referred  as  hyperbolic  balance  laws)  has  attracted  much  attention. 
In  one  space  dimension,  balance  laws  take  the  form  of 

Ut  +  f{U),  =  s{U,x),  (1.1) 

where  U  is  the  solution  vector  with  the  corresponding  flux  f{U),  and  s{U,x)  is  the  source 
term.  They  usually  admit  steady  state  solutions,  in  which  the  source  term  is  exactly  bal¬ 
anced  by  the  flux  gradient.  A  challenge  in  the  numerical  analysis  of  such  balance  laws  is 
to  maintain  these  steady  states,  and  to  compute  their  perturbations  accurately.  Indeed,  if 
a  scheme  cannot  balance  the  effects  of  convective  fluxes  and  source  term,  it  may  introduce 
spurious  oscillations  near  equilibria,  unless  the  mesh  size  is  extremely  rehned.  To  save  the 
computational  cost,  well-balanced  methods,  which  preserve  exactly  these  steady-state  so¬ 
lutions  up  to  machine  accuracy,  are  specially  designed  and  work  well  on  relatively  coarse 
meshes.  A  typical  example  considered  extensively  in  the  literature  for  balance  laws  is  the 
shallow  water  equation  with  a  non-flat  bottom  topology.  Many  researchers  have  developed 
well-balanced  methods  for  the  shallow  water  equation  using  different  approaches,  see,  e.g. 
[4,  1,  8,  2,  13,  22,  20]  and  the  references  therein. 

Another  important  example  of  the  hyperbolic  balance  laws  is  the  gas  dynamics  system  un¬ 
der  gravitational  held.  Near  the  equilibrium  state,  there  exists  the  hydrostatic  balance  where 
the  hux  produced  by  the  pressure  is  canceled  by  the  gravitational  source  term.  Many  astro- 
physical  problems  involve  the  hydro  dynamical  evolution  in  a  gravitational  held,  therefore  it 
is  essential  to  correctly  capture  the  ehect  of  gravitational  force  in  the  simulations,  especially 
if  a  long-time  integration  is  involved,  for  example  in  modeling  star  and  galaxy  formation. 
Improper  treatment  of  the  gravitational  force  can  lead  to  a  solution  which  either  oscillates 
around  the  equilibrium,  or  deviates  from  the  equilibrium  after  a  long  time  run.  There  have 
been  several  attempts  in  designing  well-balanced  methods  for  the  gas  dynamics,  which  take 
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care  of  the  implementation  of  the  gravitational  held.  LeVeque  and  Bale  [9]  proposed  the 
quasi-steady  wave-propagation  methods  for  an  ideal  gas  subject  to  a  static  gravitational 
held.  A  Riemann  problem  is  introduced  in  the  center  of  each  grid  cell  such  that  the  hux 
diherence  exactly  cancels  the  source  term.  Zingale  et  ah  [6]  investigated  the  process  of  map¬ 
ping  an  astrophysical  initial  model  from  a  stellar  evolution  code  onto  the  computational  grid 
of  an  explicit  code  while  maintaining  hydrostatic  equilibrium.  A  diherent  strategy  for  the 
construction  of  well-balanced  discretizations  with  respect  to  dominant  hydrostatics  has  been 
proposed  by  Botta  et  ah  [5]  for  the  nearly  hydrostatic  hows  belonging  to  a  certain  class  of 
solutions.  More  recently,  Xu  and  his  collaborators  [19,  23,  12]  have  extended  the  gas-kinetic 
scheme  to  the  multidimensional  gas  dynamic  equations  to  develop  well-balanced  numerical 
methods,  where  the  gravitational  potential  was  modeled  as  a  piecewise  step  function  with  a 
potential  jump  at  the  cell  interface. 

Most  of  the  works  mentioned  above  are  for  numerical  schemes  of  at  most  second  order 
spatial  accuracy.  The  main  objective  of  this  paper  is  to  design  a  hnite  diherence  WENO 
scheme  for  the  gas  dynamic  equations  with  gravitational  source  terms,  which  maintains  the 
well-balanced  property  for  the  hydrostatic  balance,  and  at  the  same  time  is  genuinely  high 
order  accurate  for  the  general  solutions.  We  introduce  a  diherent  treatment  of  the  source 
term,  which  mimics  the  WENO  approximation  to  the  hux  term,  so  that  the  exact  balance 
between  the  source  term  and  the  hux  can  be  achieved  at  the  steady  state.  The  proposed 
method  is  a  generalization  of  the  approach  to  develop  well-balanced  methods  for  the  shallow 
water  equations  and  other  balanced  laws  in  [20,  21].  The  specihc  WENO  scheme  we  use  is 
the  hfth  order  hnite  diherence  scheme  introduced  by  Jiang  and  Shu  [7].  It  uses  a  convex 
combination  of  three  candidate  stencils,  each  producing  a  third  order  accurate  hux,  to  obtain 
hfth  order  accuracy  and  an  essentially  nonoscillatory  shock  transition.  Time  discretization 
can  be  implemented  by  the  TVD  Range-Kutta  method  [17].  WENO  schemes  were  hrst 
introduced  by  Liu,  Osher  and  Chan  in  [11].  For  more  details  of  WENO  schemes,  we  refer  to 
[7,  15,  16]. 
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The  outline  of  the  paper  is  as  follows:  In  Section  2,  the  model  and  its  steady  state 
solutions  are  described.  In  Section  3,  we  develop  the  well-balanced  hnite  difference  WENO 
scheme,  which  at  the  same  time  is  genuinely  high  order  accurate  for  the  general  solutions 
of  the  gas  dynamics  equations.  Extension  to  two  dimensional  problems  is  introduced  in 
Section  4.  In  Section  5,  we  show  selective  numerical  results  in  one  and  two  dimensions  to 
demonstrate  the  behavior  of  the  proposed  hnite  difference  WENO  methods,  verifying  high 
order  accuracy,  the  well-balanced  property,  and  good  resolution  for  smooth  and  discontinuous 
solutions.  Concluding  remarks  are  given  in  Section  6. 

2  The  model 

We  consider  the  equations  governing  the  conservation  of  mass,  momentum  and  energy  of  an 
inviscid,  non-heat  conducting,  isotropic  huid.  These  gas  dynamic  equations,  coupled  with  a 
static  gravitational  potential,  are  given  by 

Pt  +  {pu)x  =  0 

{pu)t  +  {pu^  +p)x  =  -p(t>x 

Et  + {{E +  p)u)a,  = -pucpa^,  (2.1) 

in  one  space  dimension,  where  p  denotes  the  huid  density,  u  is  the  velocity,  p  represents  the 
pressure,  and  E  =  \pu^  +p/ (py  —  1)  is  the  non-gravitational  energy  which  includes  the  kinetic 
and  internal  energy  of  the  huid.  7  is  the  ratio  of  specihc  heats  and  0  =  0(x)  is  the  time 
independent  gravitational  potential. 

2.1  Hydrostatic  balance 

For  the  static  gravitational  potential  0(x),  we  are  interested  in  preserving  the  hydrostatic 
balance,  which  is  a  special  steady  state  solution  to  (2.1) 

p  =  p{x),  n  =  0,  p^  =  -p(t)x,  (2.2) 
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with  a  constant  temperature  and  zero  velocity.  For  an  ideal  gas,  we  have  the  relation 

P  =  pRT,  (2.3) 

where  R  is  the  gas  constant,  and  T  is  the  temperature.  Combined  with  the  steady  state 
equation  px  =  —p4>x  from  (2.2),  it  becomes 

p{x)  =  poexp  ,  (2.4) 

which  leads  to  the  special  steady  state 

p  =  poexp  ,  u  =  0,  p  =  RTp  =  RTpoexp(^-^^  ,  (2.5) 

with  constant  temperature  T. 

The  simplest  and  most  commonly  encountered  case  is  the  linear  gravitational  potential 
held:  d(f)/dx  =  g,  with  the  corresponding  hydrostatic  balance 

p  =  Poexp{-gpox/po),  u  =  0,  p  =  poexp{-gpox/po).  (2.6) 

3  Well-balanced  finite  difference  WENO  methods 

In  this  section  we  design  a  genuinely  high  order  accurate  well-balanced  WENO  method 
for  the  gas  dynamic  equations  (2.1)  with  gravitational  source  terms.  The  key  idea  is  to 
discretize  the  source  term  by  a  hnite  difference  WENO  formula  consistent  with  that  for  the 
hux.  We  will  concentrate  our  discussion  on  the  one-dimensional  case.  Generalization  to  the 
two  dimensional  problems  will  be  presented  in  Section  4. 

Well-balanced  methods  are  specially  designed  to  preserve  exactly  the  steady-state  so¬ 
lutions  up  to  machine  accuracy  with  relatively  coarse  meshes.  High-order  well-balanced 
WENO  methods  have  been  designed  for  the  shallow- water  equations  by  the  same  authors  in 
[20].  The  main  idea  there  is  to  decompose  the  source  term  into  a  sum  of  two  terms  first,  and 
discretize  each  of  them  independently  using  a  finite  difference  formula  consistent  with  that 
of  approximating  the  flux  derivative  terms  in  the  conservation  law.  The  same  technique  has 
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been  generalized  to  other  hyperbolic  balance  laws  in  [21]  and  to  hybrid  WENO  schemes  in 
[10].  Similar  idea  has  also  been  employed  in  designing  well-balanced  method  for  the  Euler 
equation  in  [5],  where  the  same  discrete  gradient  operator  has  been  used  to  approximate 
the  pressure  gradient  and  gravitational  potential  gradient.  In  this  paper,  we  generalize  this 
idea  to  the  gas  dynamic  equations  to  design  a  genuinely  high  order  well-balanced  WENO 
method.  We  would  like  to  preserve  exactly  the  steady  state  solution  satisfying  (2.5).  The 
special  steady  state  of  the  form  (2.6) 

p  =  cexp(— gfx),  M  =  0,  p  =  cexp(— gfx),  (3.1) 

coupled  with  the  linear  gravitational  potential  field 

=  g,  (3.2) 

will  be  used  as  an  example  to  illustrate  the  basic  idea  in  this  section.  More  general  equilib¬ 
rium  (2.5)  can  be  handled  following  the  same  technique,  and  will  be  explained  briefly  at  the 
end  of  this  section. 

The  hrst  step  in  designing  the  well-balanced  method  is  to  rewrite  the  source  term  following 
the  guideline  presented  in  [21].  The  gas  dynamic  equations  (2.1)  can  be  reformulated  as 

Pt  +  {pu)x  =  0 

{pu)t  +  {pu^  +p)x  =  p  eyiv{gx){ey.p{-gx))x 

Et  -h  {{E  +  p)u)x  =  pueyip{gx){fiyip{-gx))x,  (3.3) 

where  the  gravitational  source  —pg  is  replaced  by  p eyip{gx){isyip{—gx))x-i  and  —pug  is  treated 
in  the  same  way.  The  main  motivation  of  such  a  change  is  to  let  the  source  term  and  the 
corresponding  flux  term  share  similar  form  in  the  case  of  the  steady  state  solution  (3.2). 
As  we  will  see  below,  this  new  form  of  the  source  term  is  crucial  for  the  design  of  our 
well-balanced  method.  We  refer  to  [21]  for  more  motivation  and  explanation,  as  well  as  the 
extension  to  more  general  source  terms. 
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For  simplicity,  we  denote  the  one-dimensional  Euler  equations  (3.3)  by 


Ut  +  f{U),  =  s{U,(l)), 


where  U  =  {h,  hu,  E)'^  with  the  superscript  T  denoting  the  transpose,  f{U)  is  the  flux  and 
s{U,(j))  =  {0,  —p(j)x,  pu4>x)'^  is  the  source  term.  We  also  assume  the  mesh  is  uniform  with 
mesh  size  Ax.  In  a  finite  difference  scheme,  our  computational  variables  are  Uj{t),  which 
approximate  the  point  value  at  xj.  The  finite  difference  WENO  scheme  is  given  by 


A 

dt 


U,(t)  + 


1 

Ai, 


I, ((7), 


(3.4) 


with  being  the  numerical  flux,  and  Sj{U)  being  a  high  order  approximation  to  the  source 
term  at  the  point  Xj. 

A  well-balanced  method  is  one  that  balances  the  flux  and  source  term  exactly,  i.e. 
{jj+i  — /j_i  j  /Axj  =  Sj{U),  at  the  steady  state  solution  (3.1).  Following  the  ideas  in 
[20],  we  will  start  by  considering  an  identical  linear  finite  difference  operator  for  the  flux 
derivative  and  the  derivatives  in  the  source  terms.  A  linear  finite  difference  operator  D  is 
defined  to  be  one  satisfying  D{agi  +  6(72)  =  ciD{gi)  +  bD{g2)  for  constants  a,  b  and  arbitrary 
grid  functions  gi  and  g2-  A  scheme  for  (3.3)  is  said  to  be  a  linear  scheme  if  all  the  spa¬ 
tial  derivatives  are  approximated  by  linear  finite  difference  operators.  At  the  steady  state 
solution  (3.1),  for  any  consistent  linear  scheme,  the  first  equation  {pu)x  =  0  and  the  third 
equation  ((E  +p)u)x  =  puexp{gx){exp{—gx))x  are  both  satisfied  exactly  since  the  velocity 
u  equals  to  zero.  The  truncation  error  for  the  second  equation  takes  the  form  of 


Diipu^  +p)  -  pexp(gx)D2(exp(-gx)), 


(3.5) 


which  reduces  to 


Di(p)  -  pexp(gx)D2(exp(-gx)), 

where  Ei  and  D2  are  linear  finite  difference  operators  used  to  approximate  the  spatial  deriva¬ 
tives.  For  any  linear  scheme  which  also  satisfies  that  Ei  =  D2  for  the  steady  state  solution 
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(3.1),  the  truncation  error  for  the  second  equation  reduces  to 


Di{p)  -  pexp{gx)D2{exp{-gx))  =  Di{cexp{-gx))  -  cexp{-gx)  exp{gx)D2{exp{-gx))  =  0, 

at  the  steady  state  solution  (3.1),  which  presents  the  desired  well-balanced  property.  There¬ 
fore,  we  have 

Proposition  3.1:  For  the  gas  dynamic  equations  (2.1)  with  the  linear  gravitational  potential 
held  (3.2),  linear  scheme  which  satishes  Di  =  D2  for  the  steady  state  solution  (3.1)  are  well- 
balanced,  i.e.  it  can  preserve  these  steady  state  solution  exactly. 

We  now  extend  the  well-balanced  property  to  high  order  nonlinear  finite  difference  WENO 
scheme  [7,  3],  in  which  the  nonlinearity  comes  from  the  nonlinear  weights  and  the  smooth 
indicators.  We  would  like  to  make  minor  modihcations  to  the  WENO  scheme,  so  that  the 
well-balanced  property  is  maintained  without  affecting  the  accuracy  and  nonlinear  stability. 

To  simplify  the  presentation,  we  hrst  consider  the  hnite  difference  WENO  scheme  without 
a  flux  splitting  (e.g.  the  WENO-Roe  scheme  as  described  in  [7])  or  the  local  characteristic 
decomposition.  At  the  steady  state  solution  (3.1),  the  first  and  third  equations  in  (3.3)  can 
be  satished  exactly,  since  u  =  0  and  the  WENO  approximation  to  {pu)^  and  {{E  +  p)u)x  is 
consistent.  There  are  two  derivatives  in  the  second  momentum  equation  of  (3.3):  {pu^  +p)x 
and  (exp(— ^fa:))^;.  The  WENO  approximation  to  the  flux  derivative  term  {pu^+p)x  proceeds 
as  usual.  Notice  that  the  WENO  approximation  to  where  d  =  pv?  +p  can  be  eventually 
written  out  as 

r 

^^\x=xj  ~  ^  ^  O'kdk+j  =  Dd{d)j  (3-6) 

k=—r 

where  r  =  3  for  the  fifth  order  WENO  approximation  and  the  coefficients  depend  non- 
linearly  on  the  smoothness  indicators  involving  the  grid  function  d.  For  the  other  derivative 
(exp(— ^fx))^;  in  the  source  term,  as  explained  in  [21],  the  key  idea  now  is  to  use  the  finite 
difference  operator  with  d  =  pu^+p  fixed,  namely  to  use  the  same  coefficients  obtained 
through  the  smoothness  indicators  of  d  in  (3.6),  and  apply  it  to  approximate  (exp(— ^fx))^;. 


Thus 


~  flfc  exp(-^Xfc+j)  =  Dd{ew{-9x))j. 

k=—r 

Clearly,  with  d  =  f{u,  x)  being  fixed,  the  finite  difference  operator  defined  from  the  high 
order  WENO  procedure,  is  a  high  order  accurate  /mear  approximation  to  the  first  derivative 
for  any  grid  function.  Therefore,  the  overall  high  order  accuracy  will  not  be  affected,  and 
the  proof  for  Propositions  3.1  will  go  through.  We  can  conclude  that  the  high  order  finite 
difference  WENO  scheme  as  stated  above,  without  the  flux  splitting  or  local  characteristic 
decomposition,  and  with  the  special  handling  of  the  source  terms  described  above,  maintains 
exactly  the  steady  state  solution  (3.1). 

Now,  we  consider  WENO  scheme  with  the  local  characteristic  decomposition,  which 
is  typically  used  in  high  order  WENO  scheme  for  system  to  obtain  better  non-oscillatory 
properties  for  strong  discontinuities.  To  compute  the  numerical  flux  at  we  first  compute 
an  average  state  between  Ui  and  Uj+i,  using  either  the  simple  arithmetic  mean  or  a  Roe’s 
average  [14].  The  local  characteristic  matrix  R,  consisting  of  the  right  eigenvectors  of  the 
Jacobian  at  Wj+i,  is  then  fixed.  The  neighboring  point  values  of  the  flux  functions  needed 
for  computing  the  numerical  flux  is  projected  to  the  local  characteristic  fields  determined  by 
R~^.  The  numerical  fluxes  are  computed  in  the  characteristic  direction,  and  then  projected 
back  into  the  physical  space  by  left  multiplying  with  R,  yielding  finally  the  numerical  fluxes 
in  the  physical  space.  In  this  process,  we  notice  that  (3.6)  still  holds,  while  d  now  becomes  a 
vector  grid  function  (pn,  pn^  +  p,  (E  +  p)u)'^ ,  and  are  3x3  matrices  depending  nonlinearly 
on  the  smoothness  indicators  involving  the  grid  function  d.  The  key  idea  is  still  to  use  the 
difference  operator  with  d  =  (pu,  pu^  +  p,  (E  +  p)u)'^  fixed,  and  apply  it  to  approximate 
(0, exp(— px), exp(— px))^  in  the  source  terms.  The  remaining  arguments  stay  the  same  as 
above,  and  we  can  prove  that  after  the  incorporation  of  a  local  characteristic  decomposition, 
the  WENO  scheme  still  maintains  exactly  the  steady  state  solution. 

At  the  end,  WENO  scheme  with  a  Lax-Friedrichs  flux  splitting,  such  as  the  WENO-LF 
and  WENO-LLF  scheme  described  in  [7],  is  considered.  Here  the  flux  f{U)  is  split  into 
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f~^{U)  and  /  (U),  defined  by 


fHu) 


pu  \  (  P  \ 

pv?  +p  ±  Oj  pn  , 

{E+p)ul  \Ej_ 


(3.7) 


for  the  i-th  characteristic  held,  where  a*  =  niax„  |Ai(M)|  with  Xi{u)  being  the  zth  eigenvalue 
of  the  Jacobian  f\U),  and  the  maximum  is  taken  over  either  a  local  region  (for  WENO-LLF 
method)  or  a  global  region  (WENO-LF).  The  artihcial  viscosity  terms  Eai{p,  pu,  E)'^  are 
important  to  provide  a  non-oscillatory  solution,  however  they  will  destroy  the  well-balanced 


property.  We  now  modify  this  hux  splitting  to 


fHu) 


pu  \  /  pexp{gx)  \ 

pu^+p  ±  a'  puexp{gx)  , 
{E  +  p)u  j  y  Eexp{gx)  j 


(3.8) 


by  replacing  ±Q;j(p,  pu,  E)^  with  ±Q!'(pexp(px),  pMexp(pa:),  E'exp(pa:))^.  The  coefficient  a' 


is  given  by 


a[  =  CKj  maxexp(— px), 

X 


(3,9) 


where  the  maximum  is  taken  over  a  local  or  global  region,  in  order  to  maintain  enough 
artihcial  viscosity.  This  modihcation  does  not  ahect  accuracy,  which  relies  only  on  the  fact 
f{U)  =  f^{U)  +  f~{U).  Our  motivation  of  using  ±Q;'(pexp(pa;),  pn  exp(px),  Eexp(px))^ 
comes  from  the  fact  that  they  become  constant  vectors  in  the  case  of  steady  state  solution 
(3.1).  Thus,  by  the  consistency  of  the  WENO  approximation,  their  contribution  towards 
the  numerical  hux  approximation  is  zero.  The  hux  splitting  WENO  approximation  in  this 
situation  becomes  simply  f^{U)  =  ^f{U),  hence  the  steady  state  solution  is  preserved  as 
before,  if  we  simply  split  the  derivatives  in  the  source  term  as: 


(  ° 

exp(-pa;) 

V  exp(-pa:) 


0 

exp(-pa;) 

exp(-pa;) 


0 

exp(-px) 

exp(-px) 


(3.10) 


and  apply  the  same  hux  splitting  WENO  procedure  to  approximate  them  with  the  nonlinear 


coefficients  coming  from  the  WENO  approximations  to  f^{U)  respectively.  We  have  thus 


proved  that 


10 


Proposition  3.2:  The  WENO-Roe,  WENO-LF  and  WENO-LLF  schemes  as  stated  above 
are  well-balanced  for  the  steady  state  solution  (3.1)  and  can  maintain  the  original  high  order 
accuracy.  □ 


At  the  end,  we  present  how  the  high  order  well-balanced  WENO  method  can  be  designed 
for  the  more  general  steady  state  (2.5).  We  hrst  rewrite  the  gas  dynamic  equations  (2.1)  as 


Pt  +  {pu)x  =  0 

{pu)t  +  {pv^  +v)x  =  RTpexp{^)  (^exp(-^) 
Et  +  {{E  +  p)u)x  =  RTpuexp{^)  (exp{-^] 


(3.11) 


Following  exactly  the  same  technique  as  stated  above,  with  the  flux  splitting  (3.8)  replaced 
by 


fHu)  = 


pu 


pexp( 


RT) 


pu^+p  io'  pwexp(;^) 
{E+p)u  )  \  Eexp{^) 


(3.12) 


and  the  source  term  splitting  (3.10)  replaced  by 


0 


0 


^M-4r)  =2  ) 

)  /  X  V  exp(-;j^) 


0 

T  I 

exp{-^) 


(3.13) 


we  can  show  that  the  resulting  WENO  method  is  both  high  order  accurate  and  well-balanced. 


Total  variation  diminishing  (TVD)  high  order  Runge-Kutta  time  discretization  [17]  is 
used  in  practice  for  stability  and  to  increase  temporal  accuracy.  For  example,  the  third 
order  TVD  Runge-Kutta  method  is  used  in  the  simulation  in  this  paper: 

^(1)  _  V”  +  AtJ^(V”)  (3.14) 

=  -f/"  +  - +  AtJ^(t/(2)))  , 

3  3 

where  E{U)  is  the  spatial  operator. 
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4  Two-dimensional  extension 


In  this  section,  we  construct  a  well-balanced  WENO  scheme  on  rectangular  meshes  to  solve 
the  two-dimensional  gas  dynamic  equations  with  the  gravitational  field,  which  take  the  form 


Pt  +  +  {,pv)y  =  0 

{pu)t  +  {pu^  +  p)x  +  {puv)y  =  -p(t)x 

{pv)t  {puv)x  +  {pv^  +  p)y  =  -p(t)y 

Et  +  {{E  +  p)u):x  +  {{E  -h  p)v)y  =  -pu(t)x  -  pv(j)y,  (4.1) 


where  {u,v)  is  the  velocity  of  the  fluid,  and  p,  p,  follow  the  definitions  below  (2.1). 

E  =  \p{u^  +  n^)  +p/(7  —  1)  is  the  non-gravitational  energy.  The  hydrostatic  balance  we  are 
interested  to  preserve  is  the  constant  temperature  and  zero  velocity  steady  state  solution: 


RT  )  ’ 

and  the  steady  state  solution  corresponding  to  the  linear  gravitational  potential  field:  d(f)/dx  = 
gi  and  dcfy/dy  =  g2  takes  the  form 


(4.2) 


p  =  po  exp 


RT 


M  =  n  =  0,  p  =  RT p  =  RT po  exp 


p  =  Poexp{-po{gix  +  g2y)/po),  u  =  v  =  0,  p  =  poexp{-po{gix  +  g2y)/po).  (4.3) 


It  is  straightforward  to  extend  the  high  order  finite  difference  WENO  scheme  to  multiple 
space  dimensions,  by  simply  approximating  each  spatial  derivative  along  the  relevant  coor¬ 
dinate.  This  is  one  of  its  major  advantages  over  the  finite  volume  method.  The  conservative 
approximation  to  the  derivative  from  point  values  is  as  simple  in  multi-dimensions  as  in 
one  dimension.  In  fact,  for  fixed  j,  if  we  take  W{x)  =  f{U{x,yj)),  then  we  only  need  to 
perform  the  one-dimensional  WENO  approximation  to  W(x)  to  obtain  an  approximation  to 
W\xi)  =  fx{,U{xi,  pj)).  See  again  [7,  15]  for  more  details  of  finite  difference  WENO  schemes 
in  multi-dimensions. 

Following  the  steps  in  designing  a  well-balanced  WENO  scheme  for  the  one-dimensional 
gas  dynamic  equations  in  Section  3,  we  find  out  that  it  is  also  straightforward  to  extend 
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these  results  to  two  dimensions.  For  the  steady  state  solution  (4.2),  we  rewrite  the  source 
terms  as 


in  (4.1),  and  the  one  dimension  procedure  described  in  Section  3  is  followed  in  each  of  the  x 
and  y  directions.  The  residues  are  then  summed  up  and  the  time  discretization  is  still  by  a 
TVD  Runge-Kutta  method  (3.14).  All  the  desired  properties  proved  in  the  one-dimensional 
case,  such  as  high  order  accuracy  and  the  well-balanced  property,  are  still  valid  in  the  two 
dimensional  case.  We  restricted  our  discussion  in  this  section  to  two  space  dimensions  only, 
although  the  algorithm  can  be  easily  designed  for  three  dimension  as  well. 

5  Numerical  examples 

In  this  section  we  present  numerical  results  of  our  fifth  order  well-balanced  finite  difference 
WENO  methods  for  the  one-  and  two-dimensional  gas  dynamic  equations  with  gravitational 
source  terms.  Time  discretization  is  by  the  third  order  TVD  Runge-Kutta  time  discretization 
(3.14).  Unless  otherwise  specihed,  the  CFL  number  is  taken  as  0.6. 

5.1  Shock  tube  under  gravitational  field 

The  first  test  case  is  the  standard  Sod  test,  coupled  with  the  gravitational  held.  Following 
the  problem  setup  in  [12],  the  computational  domain  is  set  as  [0,  1],  and  the  initial  conditions 
are  given  by 


p  =  1,  M  =  0,  p  =  1,  if  a;  <  0.5, 


p  =  0.125,  M  =  0,  p  =  0.1,  ifa;>0.5. 
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The  gravitational  field  0  takes  a  value  oi  g  =  (px  =  1,  and  7  =  1.4.  We  compute  this  problem 
using  our  well-balanced  finite  difference  WENO  method  with  reffection  boundary  conditions 
and  100  uniform  meshes.  The  solutions  at  time  t  =  0.2  are  shown  in  Figure  5.1  for  the 
density,  velocity,  pressure  and  the  energy.  We  also  plot  a  reference  solution,  computed  by 
the  traditional  finite  difference  WENO  method  with  much  refined  2000  uniform  meshes,  in 
these  figures  to  provide  a  comparison.  Due  to  the  gravitational  force,  the  density  distribution 
is  pulling  towards  the  left  direction,  and  negative  velocity  appears  in  some  regions.  By 
comparing  the  results  in  these  figures,  we  observe  that  this  problem  is  well  solved  by  the 
proposed  numerical  method  on  the  relatively  coarse  mesh  of  100  mesh  points. 

5.2  One  dimensional  isothermal  equilibrium  solution 

This  test  case  was  first  proposed  by  LeVeque  and  Bale  [9]  and  later  used  in  [19,  12],  to 
demonstrate  the  capability  of  the  proposed  scheme  for  computations  on  the  small  pertur¬ 
bation  of  a  steady  state.  The  computational  domain  is  set  as  [0,1].  We  consider  an  ideal 
gas  with  7  =  1.4  and  the  linear  gravitational  field  (px  =  9,  which  stays  at  an  isothermal 
equilibrium  solution  in  the  form  of  (2.6), 

Po{x)  =  po{x)  =  exp(— x),  and  Uo{x)  =  0.  (5.1) 

The  gravitational  force,  with  g  =  (px  =  acts  in  the  negative  x  direction. 

5.2.1  Well-balanced  test 

We  first  show  an  example  to  demonstrate  the  well-balanced  property  of  the  scheme.  The 
initial  condition  is  taken  as  the  steady  state  solution  (5.1)  which  should  be  exactly  preserved. 
We  compute  the  solution  until  t  =  2  using  200  uniform  mesh  points.  In  order  to  demonstrate 
that  the  steady  state  is  indeed  maintained  up  to  round-off  error,  we  use  single  precision  and 
double  precision  to  perform  the  computation,  and  show  the  errors  of  p,  pu  and  E  in 
Table  5.1.  We  can  clearly  see  that  the  errors  are  at  the  level  of  round-off  errors  for  different 
precisions,  which  verifies  the  well-balanced  property. 
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Figure  5.1:  The  numerical  solutions  of  the  Sod  test  under  gravitational  held  in  Section  5.1 
at  time  t  =  0.2.  Top  left:  density  distribution;  Top  right:  velocity  distribution;  Bottom  left: 
energy  distribution;  Bottom  right:  pressure  distribution. 
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Table  5.1;  error  for  different  precisions  for  the  steady  state  solution  (5.1)  in  Section  5.2. 


precision 

p  pu  E 

single 

2.65E-07  8.32E-08  1.37E-07 

double 

1.15E-15  1.43E-16  5.05E-16 

5.2.2  Perturbation  of  the  equilibrium  solution 

We  now  impose  a  small  perturbation  to  the  initial  pressure  state 

p{x,  f  =  0)  =  po{x)  +  rj  exp(— 100(a;  —  0.5)^),  (5.2) 

where  r;  is  a  non-zero  perturbation  constant.  Two  cases  have  been  run:  rj  =  0.01  and 
7]  =  0.0001.  For  small  p,  this  disturbance  should  split  into  two  waves,  propagating  left  and 
right.  We  compute  the  solution  until  t  =  0.25  with  200  grid  points  and  simple  transmissive 
boundary  conditions.  The  results  are  shown  in  Figure  5.2.  The  initial  pressure  perturbation 
is  included  as  the  dashed  line.  For  the  purpose  of  comparison,  we  also  plot  a  reference 
solution,  obtained  with  2000  points.  The  traditional  non-well-balanced  WENO  method, 
with  the  straightforward  calculation  of  the  source  terms,  has  been  tested  for  these  two  cases, 
and  their  results  are  also  included  in  Figure  5.2.  We  notice  that,  the  traditional  WENO 
methods  are  able  to  capture  the  big  perturbation  well,  but  do  not  perform  well  for  the 
small  perturbation  test  case  with  a  coarse  mesh  of  200  mesh  points.  With  the  well-balanced 
technique,  these  small  perturbations  are  well  captured  on  this  coarse  mesh. 

5.3  One  dimensional  gas  falling  into  a  fixed  external  potential 

We  consider  the  same  gas  dynamic  equations  with  static  gravitational  potential  (2.1)  and 
a  constant  temperature  and  zero  velocity  steady  state  solution  (2.5).  The  following  steady 
state  is  taken  from  the  paper  by  Slyz  and  Prendergast  [18],  and  has  also  been  considered  in 
[19,  12]. 
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Figure  5.2:  The  pressure  perturbation  of  a  hydrostatic  solution.  The  result  of  the  well- 
balanced  WENO  method  with  200  and  2000  grid  points,  and  that  of  the  non-well-balanced 
(denoted  by  non-wb)  WENO  method  with  200  grid  points.  Left:  rj  =  0.01  in  (5.2);  Right: 

T]  =  0.0001. 

The  gravitational  potential  has  the  form  of  a  sine  wave, 

0(a;)  =  sin  (5.3) 

where  L  is  the  computational  domain  length  and  00  is  the  amplitude.  The  steady  state  takes 
the  form  of 


f  0 

P  =  Poexp(^-  — 
with  a  constant  temperature  T. 

5.3.1  Well-balanced  test 


u  =  0, 


p  =  RT pQ  exp 


(5.4) 


We  start  with  a  test  to  verify  the  well-balanced  property.  Consider  an  ideal  gas  with  7  =  5/3, 
and  the  parameters  po  =  1,  R  =  1,  T  =  0.6866,  L  =  64,  0o  =  0.02  in  (5.3)  -  (5.4).  The  initial 
condition  is  taken  as  the  steady  state  (5.4),  and  we  compute  the  solution  until  t  =  50  using 
200  uniform  mesh  points.  In  order  to  demonstrate  that  the  steady  state  is  maintained  up 
to  round-off  error,  we  use  single  precision  and  double  precision  to  perform  the  computation, 
and  show  the  errors  of  p,  pu  and  E  in  Table  5.2.  We  can  easily  observe  that  the  errors 
are  at  the  level  of  round-off  errors  for  different  precisions. 
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Table  5.2:  error  for  different  precisions  for  the  steady  state  solution  (5.4)  in  Section  5.3. 


precision 

p  pu  E 

single 

1.52E-07  9.63E-08  1.79E-07 

double 

1.24E-15  1.93E-16  3.94E-16 

5.3.2  Convergence  test 


In  this  test,  we  start  from  an  initial  state,  which  is  a  small  perturbation  of  the  steady  state 
(5.4),  and  will  eventually  converge  to  an  isothermal  hydrostatic  distribution  after  long  time. 
The  initial  condition  is  given  by 


p  =  po  exp 


u  =  0, 


p  =  RT pq  exp 


+  0.001  exp(-10(x 


32)"), 

(5.5) 


where  the  parameters  i?,  T  et  ah  are  given  in  Subsection  5.3.1.  We  run  the  simulation 
with  64  uniform  mesh  points  for  1,000,000  time  steps,  and  plot  the  numerical  results  at 
the  hnal  time  in  Figure  5.3.  As  a  comparison,  we  also  include  the  numerical  results  of  the 
traditional  non- well-balanced  WENO  method,  where  the  constant  temperature  distribution 
of  the  steady  state  solution  is  not  well-captured  on  this  coarse  mesh. 


5.4  Testing  the  orders  of  accuracy 

In  this  example  we  check  the  numerical  orders  of  accuracy  when  the  schemes  are  applied  to 
the  following  two  dimensional  problem.  Consider  the  equations  (4.1)  with  a  linear  gravita¬ 
tional  held  0a;  =  =  1.  An  exact  solution  takes  the  form  of 

p{x,  yR)  =  l  +  0.2  sin(7r(a;  +  y-  t{uo  +  Vq))), 
u{x,y,t)  =  Uo,  v{x,y,t)  =  Vo, 

p{x,  y,  t)  =po  +  t{uo  +  Vo)  -  X  -  y  +  0.2  cos(7r(a;  +  y-  t{uo  +  Vo)))/7r, 

in  the  domain  [0,  2]  x  [0,  2].  uq  =  vo  =  1  and  po  =  4.5  are  chosen  in  this  test.  The  CFL 
condition  is  taken  as  0.2  and  the  time  step  At  is  taken  to  be  proportional  to  (1/Ax  -|- 
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Figure  5.3:  The  numerical  solutions  of  well-balanced  method  (solid  line)  and  non-well- 
balanced  method  (square  box,  denoted  by  non-wb)  for  the  convergence  test  in  Section  5.3.2 
after  1,000,000  time  steps.  Top  left:  density  distribution;  Top  right:  velocity  distribution; 
Bottom  left:  pressure  distribution;  Bottom  right:  temperature  distribution. 
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so  that  the  temporal  accuracy  is  of  the  same  order  as  the  spatial  accuracy.  The 
exact  solutions  are  used  as  the  boundary  condition.  We  run  the  simulation  until  the  stop 
time  t  =  0.1  and  compute  the  numerical  error.  Table  5.3  contains  the  errors  and  orders 
of  accuracy.  We  can  clearly  see  that,  in  this  two  dimensional  test  case,  the  designed  high 
order  accuracy  is  achieved  for  the  hnite  difference  WENO  scheme. 


Table  5.3:  errors  and  numerical  orders  of  accuracy  for  the  example  in  Section  5.4. 


Number 

P 

pu 

pv 

E 

of  cells 

error 

order 

error 

order 

error 

order 

error 

order 

8x8 

3.88E-03 

3.66E-03 

3.66E-03 

4.83E-03 

16  X  16 

5.06E-04 

2.93 

4.85E-04 

2.91 

4.85E-04 

2.91 

5.46E-04 

3.14 

32  X  32 

2.75E-05 

4.20 

2.72E-05 

4.15 

2.72E-05 

4.15 

3.02E-05 

4.17 

64  X  64 

1.15E-06 

4.57 

1.15E-06 

4.56 

1.15E-06 

4.56 

1.30E-06 

4.53 

128  X  128 

4.33E-08 

4.74 

4.29E-08 

4.73 

4.29E-08 

4.73 

5.10E-08 

4.68 

256  X  256 

1.54E-09 

4.81 

1.57E-09 

4.78 

1.57E-09 

4.78 

1.81E-09 

4.81 

5.5  Two  dimensional  isothermal  equilibrium  solution 

The  following  test  is  used  to  demonstrate  the  well-balanced  property  and  the  capability  of 
the  proposed  methods  for  capturing  the  small  perturbation  of  an  isothermal  equilibrium 
solution  in  the  two  dimensional  case.  We  set  the  computational  domain  as  the  unit  square. 
Consider  an  ideal  gas  with  7  =  1.4  and  the  linear  gravitational  field  (j)x  =  (t>y  =  9-  The 
isothermal  equilibrium  state  under  consideration  takes  the  form  of 

P{x,y)  =  Poexp  [-  —  {x  +  y)  ]  , 
p{x,  y)  =  Po  exp  (  -  —  {x  p)  )  , 

V  Po  J 

with  the  parameters  po  =  1.21,  po  =  I  and  g  =  1- 

5.5.1  Well-balanced  test 

The  two  dimensional  well-balanced  property  is  hrst  tested.  We  take  the  initial  condition  as 
the  steady  state  solution  (5.6),  and  compute  the  solution  until  t  =  1  using  50  x  50  uniform 


u{x,y)  =  v{x,y)  =  0, 


(5.6) 
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grid  points.  In  order  to  demonstrate  that  the  steady  state  is  indeed  maintained  up  to  round¬ 
off  error,  we  use  single  precision  and  double  precision  to  perform  the  computation,  and  show 
the  errors  of  p,  pu,  pv  and  E  in  Table  5.4.  We  can  clearly  see  that  the  errors  are  at  the 
level  of  round-off  errors  for  different  precisions,  which  verifies  the  well-balanced  property. 

Table  5.4:  error  for  different  precisions  for  the  steady  state  solution  (5.6)  in  Section  5.5.1. 


precision 

P 

pu 

pv 

E 

single 

1.72E-08 

3.31E-08 

3.38E-08 

3.42E-08 

double 

1.15E-16 

9.46E-17 

9.46E-17 

1.38E-16 

5.5.2  Perturbation  of  the  equilibrium  solution 

Next,  we  study  the  capability  of  the  proposed  method  for  the  small  perturbation  to  the 
stationary  state.  Consider  a  small  perturbation  to  the  initial  pressure  state  of  (5.6), 

p(x,|/,t  =  0)  =poexp  (-  —  {x  +  y^  +17  exp  (-  {{x  -  0.3)^  +  {y  -  0-3)^))  ,  (5.7) 

\  Po  J  \  Po  ) 

centered  at  (0.3,  0.3),  where  p  is  a.  non- zero  perturbation  constant,  and  set  as  0.001  in  this 
test. 

We  compute  the  solution  until  t  =  0.15  with  50  x  50  grid  points  and  simple  transmissive 
boundary  conditions.  The  contour  plot  of  the  pressure  perturbation  is  shown  in  the  left 
plots  of  Figures  5.4  and  5.5.  We  also  run  this  test  with  the  traditional  non-well-balanced 
WENO  method,  with  the  straightforward  calculation  of  the  source  terms.  Its  result  is  shown 
in  the  right  plots  of  Figures  5.4  and  5.5.  The  density  perturbations  are  also  shown  in  Figure 
5.6.  We  notice  that,  the  non-well-balanced  WENO  method  is  not  capable  of  capturing  such 
small  perturbation  on  the  coarse  mesh,  while  the  well-balanced  one  can  resolve  it  very  well. 
We  also  run  both  methods  with  a  refined  200  x  200  mesh,  and  the  results  are  shown  in 
Figure  5.7.  As  expected,  the  results  of  non-well-balanced  WENO  schemes  improve,  and 
the  difference  between  well-balanced  and  non-well-balanced  methods  becomes  smaller  as  the 
mesh  is  rehned. 
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Figure  5.4:  The  contours  of  the  pressure  perturbation  of  a  two  dimensional  hydrostatic 
solution  in  Section  5.5.2  at  time  t  =  0.15  with  50  x  50  grid  points.  20  uniformly  spaced 
contour  lines  from  —0.0003  to  0.0003.  Left:  results  based  on  well-balanced  scheme.  Right: 
results  based  on  non-well-balanced  scheme. 


Figure  5.5:  The  3D  hgure  of  the  pressure  perturbation  of  a  two  dimensional  hydrostatic 
solution  in  Section  5.5.2  at  time  t  =  0.15  with  50  x  50  grid  points.  Left:  results  based  on 
well-balanced  scheme.  Right:  results  based  on  non-well-balanced  scheme. 
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Figure  5.6:  The  contours  of  the  density  perturbation  of  a  two  dimensional  hydrostatic  solu¬ 
tion  in  Section  5.5.2  at  time  t  =  0.15  with  50  x  50  grid  points.  20  uniformly  spaced  contour 
lines  from  —0.001  to  0.0002.  Left:  results  based  on  well-balanced  scheme.  Right:  results 
based  on  non-well-balanced  scheme. 


Figure  5.7:  The  contours  of  the  pressure  perturbation  of  a  two  dimensional  hydrostatic 
solution  in  Section  5.5.2  at  time  t  =  0.15  with  200  x  200  grid  points.  20  uniformly  spaced 
contour  lines  from  —0.0003  to  0.0003.  Left:  results  based  on  well-balanced  scheme.  Right: 
results  based  on  non-well-balanced  scheme. 
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6  Concluding  remarks 


In  this  paper  we  have  constructed  a  well-balanced  hnite  difference  WENO  scheme  of  arbitrary 
order  of  accuracy  for  the  gas  dynamics  equations  with  a  static  gravitational  held.  Special 
source  term  discretization  is  introduced  such  that  the  resulting  WENO  scheme  balances  the 
zero-velocity  and  constant  temperature  steady  state  solutions  to  machine  accuracy,  and  at 
the  same  time  maintains  the  original  high  order  accuracy  and  essentially  nonoscillatory  prop¬ 
erty  for  general  solutions.  Numerical  examples  are  given  to  demonstrate  the  well-balanced 
property,  accuracy,  good  capturing  of  the  small  perturbation  to  the  steady-state  solutions, 
and  non-oscillatory  shock  resolution  of  the  proposed  numerical  method.  Ongoing  work  in¬ 
cludes  constructing  well-balanced  methods  for  more  general  steady  state  solutions. 
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